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CONSTANT-Q  SPECTRAL  ANALYSIS  BY  MEANS  OF 
AN  APPROXIMATE  FAST  FOURIER  TRANSFORM  TECHNIQUE 

INTRODUCTION 


Constant-Q  spectral  analysis  is  an  important  signal  processing  tech- 
nique in  radar,  sonar,  speech  processing,  and  other  fields  where  spectral 
analysis  of  a signal  or  an  event  is  desired.  There  are  a plethora  of 
time-  and  frequency-domain  algorithms  to  achieve  constant-Q  spectral 
analysis.  This  report  deals  with  one  such  algorithm.  Some  of  the  theory 
behind  the  technique  is  covered  along  with  some  examples  to  give  the 
reader  a down-to-earth  understanding  of  the  technique.  A basic  knowledge 
of  spectral  analysis  and  fast  Fourier  transforms  (FFT's)  is  assumed. 

The  problem  of  implementation  of  a constant-Q  spectral  analysis 
technique,  when  real  time  spectral  analysis  is  desired,  lies  in  the 
choice  of  an  efficient  algorithm  from  the  myriad  of  available  techniques. 
Efficiency  in  this  sense  refers  to  the  joint  minimization  of  computer 
time  and  storage  while  maintaining  an  acceptable  level  of  accuracy. 

This  report  describes  an  efficient  method  of  achieving  constant-Q 
spectral  analysis.  The  algorithm  results  from  the  merging  of  two  spec- 
tral analysis  techniques: 

1.  An  approximate  FFT  technique  for  octave  processing  with  a dif- 
ferent desired  resolution  for  each  octave;  and 

2.  A constant-Q  spectral  analysis  technique  that  is  based  on 
weighted  sums  of  uniform-resolution  spectral  estimates. 

As  in  a human  marriage,  there  are  problems  in  this  merging  of  tech- 
niques because  of  subtleties  of  the  personalities  and  idiosyncracies  of 
the  two  partners.  These  problems  are  overcome  when  the  idiosyncracies 
of  the  partners  are  understood. 

This  report  is  divided  into  four  autonomous,  though  related,  sec- 
tions : 

1.  Constant-Q  octave  processing. 

2.  The  approximate  FFT  technique  for  octave  processing. 

3.  The  combination  of  the  techniques  of  1 and  2 to  achieve  con- 
stant-Q spectral  analysis  within  octaves  by  means  of  the  approximate 
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FFT  technique. 

4.  A description  and  presentation  of  initial  results  of  a FORTRAN- 
based  algorithm  on  the  Univac  1108  at  the  Naval  Underwater  Systems  Cen- 
ter (NUSC) , New  London. 

The  first  section  deals  with  constant-Q  spectral  analysis.  In  con- 
stant-Q  spectral  analysis,  the  resolution  of  spectral  estimates  is  a 
fixed  percentage  of  the  center  frequency  of  the  estimate,  and  the  center 
frequencies  of  the  spectral  estimates  are  spaced  logarithmically.  There 
are  a number  of  ways  of  achieving  constant-Q  spectral  analysis.  The 
technique  chosen  in  this  report  uses  weighted  sums  of  uniform-resolution 
spectral  estimates  to  synthesize  constant-Q  spectral  estimates  from  the 
uniform  spectTal  estimates.  The  first  section  does  not  delve  into  the 
justification  of  underlying  theory  of  the  technique;  rather,  it  covers 

1.  The  selection,  sampling,  and  quantization  of  the  weighting 
coefficients,  and 

2.  Implementation  subtleties  of  the  technique  that,  if  ignored, 
can  lead  to  serious  degradation  of  the  spectral  estimates. 

The  second  section  describes  a technique  for  performing  uniform- 
resolution  spectral  estimation  within  octaves,  based  on  an  approximate 
FFT  technique.  The  approximate  FFT  technique  is  an  efficient  technique 
that  is  especially  powerful 

1.  When  nonur.iform  spectral  resolution  is  desired  across  a fre- 
quency band  (as  in  constant-Q  spectral  analysis); 

2.  When  a zoom- like  spectral  analysis  is  desired  about  a specific 
frequency  or  group  of  frequencies;  and 

3.  When  desired  spectral  resolutions  are  so  narrow  that  the  FFT 
size  would  be  larger  than  the  memory  available. 

In  the  approximate  FFT  technique,  the  FFT  is  treated  as  a bank  of 
digital-comb  filters.  Each  filter  in  the  comb  simultaneously  provides 
the  functions  of  bandpass  filtering  and  complex  demodulation  of  the 
input  data.  Just  as  bandpass  filtering  is  normally  done  before  spectral 
analysis,  so  too  can  the  FFT  be  used  as  an  initial  bank  of  bandpass  fil- 
ters for  input  to  a second  set  of  FFT's  which  are  applied  to  each  filter 
in  the  bank.  That  is,  a short-duration  FFT  is  used  to  separate  the  in- 
put data  into  coarse  spectral  bins  and  a second  FFT  (a  vernier  FFT)  is 
applied  to  successive  outputs  of  each  coarse  spectral  bin  to  obtain  any 
desired  resolution  within  that  bin.  Since  each  bin  has  its  own  vernier, 
there  is  great  latitude  in  choosing  possible  types  of  spectral  analysis. 
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The  third  section  relates  the  techniques  of  the  first  two  sections 
and  describes  how  the  two  techniques  are  merged  to  achieve  constant-Q 
spectral  analysis  within  octaves.  The  benefits,  drawbacks,  problems, 
and  solutions  to  problems  that  are  created  by  this  merger  also  are 
covered  in  the  third  section. 

The  fourth  section  describes  the  combined  spectral -analysis  techni- 
que as  implemented  in  FORTRAN  on  the  NUSC/NL  Univac  1108.  The  program 
described  is  a specialized  algorithm,  and  suggestions  for  generalizing 
the  technique  are  given  along  with  some  initial  outputs  of  the  routine. 


CONSTANT-Q  SPECTRAL  ANALYSIS 


Constant-Q  spectral  analysis  is  called  by  many  names,  for  example, 
proportional  bandwidth  or  constant-percentage  spectral  analysis.  Because 
of  these  and  other  names,  information  under  the  heading  constant-Q  spec- 
tral analysis  is  difficult  to  find.  As  mentioned  earlier,  Constant-Q 
spectral  analysis  can  be  considered  a specialized  form  of  conventional 
spectral  analysis,  in  which  the  center  frequencies  of  the  spectral  esti- 
mates lie  along  a logarithmic  frequency  axis  and  the  bandwidth  (resolu- 
tion) of  an  estimate  is  a fixed  percentage  of  the  center  frequency  of 
that  spectral  estimate. 

Many  techniques  exist  for  implementing  constant-Q  spectral  analysis, 
and  the  technique  presented  in  this  section  is  an  efficient  method. 
Weighted  linear  combinations  of  uniform-resolution  spectral  estimates 
can  be  used  to  generate  constant-Q  spectral  estimates.  A recent  report 
by  Fred  Harris,*  of  the  Naval  Undersea  Center  (NUC) , shows  why  the 
sum  approach  works . 

Implementation  of  this  method  involves  the  following: 

1.  Breaking  an  input  data  sequence  into  octaves  and  performing 
uniform-resolution  spectral  analysis  within  each  octave;  and 

2.  Applying  a set  of  weighting  coefficients  (often  referred  to  as 
filter  coefficients)  to  the  uniform  spectral  estimates  to  obtain  con- 
stant-Q spectral  estimates. 

The  first  step  will  be  discussed  later  in  this  report,  and  this  is 
where  the  power  and  benefits  of  the  approximate  FFT  technique  will  become 
apparent.  This  section  will  discuss  the  second  step  of  this  method  and 
will  dwell  on  the  following  points: 

1.  How  does  the  weighted-sum  method  of  constant-Q  spectral  analysis 

work? 
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2.  What  octave  processing  is,  and  why  the  weighted-sum  method 
based  on  octave  processing  is  efficient. 

3.  Why  there  is  a need  for  variable  overlapping  of  the  time  sequen- 
ces prior  to  spectral  analysis. 

4.  The  selection,  digitization,  and  quantization  of  weighting 
functions  for  the  weighting  coefficients. 


{A.)  r . „ from  a length 
K|  K ' ° i \ N - 1 
|an J n = O' 


WEIGHTED  SUM  METHOD 

Assume  we  have  a set  of  spectral  estimates 

N digital  Fourier  transform  (DFT)  of  a real  input  sequence  janj  n = 0" 

Associated  with  each  A^  is  a center  frequency,  as  well  as  a bandwidth 
(or  resolution)  of  the  estimate.  Each  estimate  has  side  lobes  associated 
with  it.  These  side  lobes  are  a nuisance,  because  they  represent  inter- 
ference from  frequencies  outside  the  frequency  of  interest.  Multiplica- 
tive windowing  of  the  input  sequence,  or  convolutional  weighting  of  the 
unweighted  DFT  outputs,  is  generally  used  to  control  spurious  responses. 


Constant-Q  spectral  analysis  by  means  of  weighted  sums  represents  a 
specialized  form  of  convolutional  weighting  in  the  frequency  domain. 

This  specialized  convolutional  form  involves  the  use  of  a different  set 
of  weights  for  each  shift  in  the  convolution.  Figures  1 through  3 illu- 
strate the  principle  involved,  while  reference  1 gives  a brief  proof  of 
why  these  weighted  sums  represent  spectral  estimates  with  arbitary  cen- 
ters and  arbitrary  handwidths. 


Figures  1 and  2 display  spectral  window  functions  W(f)  and  W^(S.). 
The  function  W(f)  is  the  Fourier  transform  of  a real,  even,  temporal  - 
domain  weighting  function , and  Wj(2)  is  a sampled,  expanded,  frequency- 
shifted  version  of  W(f).  Note  that  both  the  frequency  offset  from  an 
FFT  bin  and  an  expansion  factor  can  be  specified  independently  for  each 
set  { W C JL)  > . Figure  3 displays  the  weighting  of  some  undefined  DFT  out- 
puts with  sample  spacing  A,  to  which  the  convolutional  weighting  func- 
tions are  applied.  The  center  frequency  of  the  resulting  estimate  is 
fi  = A{j.i  + 6^},  while  the  bandwidth  is  Bj  = aiB0A.  The  implications 
of  the  relationship  of  f^  and  Bj  are  important: 

1.  By  proper  choice  of  «j  and  A,  a set  of  weights  can  be  generated 
to  synthesize  proportional-bandwidth  spectral  estimates  from  uniform 
spectral  estimates  over  an  octave  of  frequencies;  and 

2.  The  same  set  of  weights  used  in  one  octave  can  be  used  for  all 
octaves  by  the  proper  choice  of  A. 
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Window  Function  Sampling,  W(f)  a Real,  Even  Window  Function 


Figure  2.  Window  Function  Sampling,  W.  (i.)  a Sampled,  Shifted,  Expanded  W(f) 
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These  two  points  account  for  the  efficiency  of  the  weighted-sum 
approach  over  other  techniques.  This  technique  saves  computer  storage, 
since  only  a single  set  of  weights  need  he  stored  for  a wide  band  of 
frequencies  because  of  the  octave-processing  characteristics  of  2,  above. 
The  computational  savings  of  the  technique  lies  in  the  choice  of  an 
efficient  algorithm  for  generating  t he  uniform  spectral  estimates  for 
each  octave. 

OVERLAPPED  PROCESSING 

A subtle,  often  overlooked,  point  of  the  weighted-sum  technique  is 
the  choice  of  overlap  of  the  input  sequences  used  to  create  the  uniform 
spectral  estimates.  With  this  constant -Q  technique,  50  percent-over- 
lapped  time  segments  are  sufficient  to  avoid  gaps  in  the  time  domain  at 
the  low  end  of  an  octave.  However,  at  the  upper  end  of  an  octave  a 75 
percent  overlap  is  needed. 

Figure  4 shows  why  this  is  true  when  contiguous  blocks  of  data  are 
processed.  Each  filter  function  has  an  effective  time  duration  that  is 
less  than  the  length  of  the  data  record  being  processed.  Figure  4a 
illustrates  the  effective  time-domain  weightings  for  a constant-Q  filter 
at  the  low  end  of  the  octave.  Note  that  for  contiguous  blocks  of  data 
the  effective  weighting  creates  gaps  in  the  data.  The  gaps  can  be 
covered  by  50  percent-overlapped  processing  (dashed-weighting  functions). 

Figure  4b  shows  the  effective  time-domain  weightings  for  a filter 
at  the  upper  end  of  an  octave.  Note  that  the  time  domain  gaps  are 
larger  and  that  a 50  percent  overlap  is  insufficient  to  cover  the  gaps. 

Theoretically,  each  proportional-bandwidth  filter  would  require  a 
different  amount  of  overlap  (ranging  between  50  percent  and  75  percent 
overlap),  depending  upon  its  location  within  the  octave  of  interest. 

This  need  for  variable  overlap  complicate  he  application  of  the  filter 
weights  and  breaks  up  the  data  rates  withp  ’verv  octave. 

WEIGHTING  FUNCTIONS 

The  generation  of  filter  coefficients  falls  into  the  domain  of 
choosing  a good  weighting  function.  The  choice  of  a good  weighting 
function  was  investigated.  The  Kai ser-Bessel  function  was  found  to  be 
convenient  and  worthwhile  and  was  used  in  this  investigation.  The 
Kai ser-Bessel  function  has  a simply-expressed  closed-form  Fourier  trans- 
form and  a main  lobe  width  and  side  lobe  levels  that  can  be  controlled 
by  a single  parameter,  B (see  figure  5). 

Appendix  A,  written  by  Dr.  A.  H.  Nuttall,  brings  out  the  following 
points  concerning  convolutional  frequency-domain  weights: 


Figure  4b.  Upper  End  of  an  Octave 
Figure  4.  Effective  Temporal  Weightings 


Figure  5.  Origin-Centered  Kaiser-Bessel  Window  Function 
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1.  Although  the  spectral  window  is,  in  general,  complex  when  the 
weighting  is  not  centered  about  the  origin,  the  spectral  window  can  be 
represented  by  real  coefficients. 

2.  Minimum  and  maximum  sampling  rates  exist  when  a limited  number 
of  samples  are  used. 

3.  The  width  of  the  main  lobe  can  be  changed  with  negligible  effect 
on  side  lobe  levels  by  varying  the  sampling  rate  within  a certain  range. 
For  seven  coefficients,  the  maximum  range  of  variation  is  2:1. 

4.  Real  filter  coefficients  can  be  used  to  obtain  constant-Q  spec- 
tral estimates  from  uniformly  resolved  spectral  estimates,  but  only  with 
the  caution  mentioned  earlier  regarding  overlapping. 


AN  APPROXIMATE  FFT  TECHNIQUE 


The  approximate  FFT  technique 

1.  Accepts  input  data  in  moderate-sized  time  segments,  as  available; 

2.  Performs  a weighted  FFT  on  each  segment  (overlapped,  if  neces- 
sary) ; 

3.  Stores  those  frequency  portions  (at  each  segment)  where  narrow- 
band  spectral  analysis  is  desired;  and 

4.  Performs  a small-size  weighted  FFT  over  the  total  data  record 
available  for  each  frequency  portion  stored. 

The  last  transform  over  time  (delay)  in  4,  above,  is  vernier  fre- 
quency analysis,  as  measured  from  the  center  of  each  bin.  It  is  not 
exact.  However,  spurious  side  lobes  can  be  controlled  adequately  by 
proper  overlapping  and  weighting  in  2 and  the  use  of  proper  phase  fac- 
tors in  4 . By  this  procedure,  we  realize  narrowband  frequency  resolution 
without  the  need  for  performing  a large-size  FFT. 

Perhaps  the  easiest  way  to  present  the  approximate  FFT  technique  is 
by  a series  of  simple  numerical  examples.  Suppose  we  segment  the  input 
data  into  1 sec  sections.  In  figure  6,  a time-frequency  diagram  illu- 
strates the  basic  operations  required.  A temporal  weighting  (e.g.,  a 
Hanning  function)  is  applied  to  each  1 sec  segment  of  input  data,  which 
is  then  subjected  to  an  FFT.  The  location  along  the  frequency  axis  of 
the  resultant  spectral  components  is  indicated  by  a vertical  line  of 
X's,  spaced  1 Hz  apart.  For  every  50  percent  shift  (1/2  sec  delay)  of 
the  temporal  weighting,  the  procedure  is  repeated  eight  times,  resulting 
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Figure  6.  Time-Frequency  Diagram 

in  a series  of  eight  complex  spectral  outputs  covering  a time  interval 
of  approximately  4 sec.  Then  the  spectral  components  in  a particular 
frequency  bin  of  interest,  as  indicated  by  the  horizontal  box  in  figure 
6,  are  (delay)  weighted  and  subjected  to  a final  8-point  (vernier)  FFT. 

For  the  example  depicting  eight  outputs  in  time,  each  separated  by 
1/2  sec,  the  resolution  of  the  vernier  FFT  is  Av  = 1/4  Hz.  Since  the 
total  frequency  coverage  of  this  vernier  FFT  is  8Av  = 2 Hz,  whereas 
the  original  separation  of  frequency  components  in  figure  6 is  only 
Af  = 1 Hz,  only  the  central  half  of  the  vernier  outputs  about  zero  fre- 
quency are  retained.  Thus,  we  can  get  complete  coverage  of  the  frequency 
scale  at  Av  = 1/4  Hz  resolution  without  having  to  conduct  the  larger 
size  FFT,  i.e.,  over  the  entire  4 sec  interval. 

To  indicate  in  detail  how  the  technique  works,  consider  a unit- 
amplitude  pure  tone  at  256  Hz  with  zero  phase,  as  shown  in  figure  7. 

To  simplify  matters,  neither  temporal  weighting  nor  overlap  are  employed. 
The  results  of  the  initial  1 sec  FFT's  are  a series  of  eight  ones  at  a 
frequency  of  256  Hz  and  zeros  at  all  other  frequencies.  Then,  when  the 
bin  centered  at  256  Hz  is  subjected  to  the  8-point  vernier  FFT,  only  the 
zeroth  location  (zero  frequency)  has  a nonzero  output.  Also,  vernier 
FFT's  conducted  on  other  bins  yield  zero  outputs  everywhere. 
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Figure  7.  Pure  Tone  at  256  Hz  (Zero  Phase)  With  Unweighted  FFT's  and  No  Overlap 
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In  figure  8,  a tone  is  placed  instead  at  256-1/2  Hz.  Now  bins  cen- 
tered at  256  and  257  Hz  have  alternating  outputs  of  magnitude  2/tt.  The 
other  bins  have  lesser  alternating  outputs  depending  on  their  proximity 
to  25b-l/2  Hz,  and  decay  according  to  a sinc-function  law  for  the 
unweighted  example  that  we  are  considering.  The  vernier  FFT  on  the  bin 
centered  at  256  Hz  now  contains  only  one  nonzero  element  at  the  fourth 
location.  A similar  condition  holds  for  other  FFT  bins,  but  with 
decreased  magnitude,  depending  on  thei r proximi ty  to  256-1/2  Hz. 

Finally,  in  figure  9 a tone  is  placed  at  256-1/8  Hz.  The  frequency 
bins  have  outputs  that  change  by  exp(i2ir/8)  for  every  1 sec  FFT,  and 
their  magnitudes  are  governed  by  a sinc-function  law.  The  vernier  FFT 
on  the  bin  at  256  Hz  now  has  only  one  nonzero  output  at  the  first  loca- 
tion. , 


We  can  see  now  how  the  proposed  technique  works.  The  location  of 
the  peak  output  of  the  vernier  FFT  for  a particular  frequency  bin  is  a 
direct  measure  of  the  tone  frequency  present  in  the  input  data,  relative 
to  the  center  of  that  particular  bin.  In  order  to  control  side  lobes  and 
spurious  responses  of  this  technique,  overlap,  temporal  weighting,  delay 
weighting,  and  phase  factors  must  be  selected  and  applied  carefully.  In 
particular,  for  50  percent  overlap  and  Dolph-Chebyshev  temporal  weight- 
ing, -33  dB  spurious  side  lobe  responses  can  he  realized.  For  75  per- 
cent overlap  and  Dolph-Chebyshev  temporal  weighting,  -86  dB  side  lobes 
are  attainable.  In  both  cases  of  overlap,  delay  weighting  is  required 
and  was  taken  as  a Hanning  function.  Although  phase  factors  were 
unnecessary  in  the  nonoverlap  example  in  figures  7 through  9,  they  are 
mandatory  when  overlap  is  used. 

For  a more  detailed  mathematical  treatment  of  how  the  technique 
works  and  for  extensions  to  other  overlaps  and  weightings,  see  reference 
2. 


CONSTANT-Q  SPECTRAL  ANALYSIS  BY  MEANS  OF  I HI  APPROXIMATE  FFT 


The  approximate  FFT  technique  can  now  be  applied  to  the  weighted- 
sum  method  for  constant-Q  spectral  analysis.  The  weighted-sum  method 
requires  uniformly  resolved,  linearly  spaced  spectral  estimates  with  a 
very  specific  format.  The  format  is  based  on  the  octave  processing 
described  earlier.  Assume  for  a moment  that  some  arbitrary  frequency 
range  has  been  broken  into  octaves  and  the  octaves  have  been  numbered, 
starting  with  the  highest  frequency  octave  as  1,  the  next  lower  octave 
as  2,  etc.  If  one  also  assumes  that  the  resolution  of  the  spectral 
estimates  within  the  upper  octave  is  v,  then  the  required  resolution  for 
any  octave  is  R(octave  #)  = y2"(°ctave  # - 1) . Another  requirement  of 
the  specific  format  is  the  variable  overlap  of  the  time  sequences. 
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Figure  10  illustrates  the  approximate  FFT  technique  implemented  as 
an  octave  processor.  This  technique  generates  with  ease  the  desired 
pattern  of  octave  resolutions.  Instead  of  selecting  single  frequencies 
from  each  temporal  FFT  (as  in  the  approximate  FFT  technique  shown  in 
figure  6),  we  now  collect  subsets  of  temporal  FFT  outputs  that  repre- 
sent octaves  of  interest.  For  example,  one  might  collect  8 of  these 
subsets  for  the  first  octave,  16  for  the  second  octave,  etc.,  and  then 
perform  a vernier  FFT  on  the  complex  output  series  from  each  coarse  FFT 
bin  within  the  subset.  Each  vernier  FFT  creates  uniformly  resolved 
spectral  estimates  that  are  then  weighted  and  summed  to  generate  con- 
stant-!} spectral  estimates.  Thus , the  proper  resolution  pattern  is 
achieved  through  doubling  the  vernier  FFT  for  each  lower  octave. 

The  approximate  FFT  technique  is  ideally  suited  to  constant -Q 
octave  processing  for  the  following  reasons  (.figure  11): 

1.  The  temporal  FFT  allows  processing  of  only  the  frequencies 
within  the  desired  octaves; 

2.  Implementation  of  the  needed  resolution  within  each  octave  is 
accomplished  by  means  of  the  second  (vernier)  FFT;  and 

3.  A limited  variable  data  overlap  can  be  implemented  by  overlap- 
ping the  temporal  FFT  outputs,  if  desired. 

There  are  indications  that  constant-Q  octave  processing  based  on 
the  approximate  FFT  technique  offers  computational  savings  over  other 
techniques . 

The  weighted-sum  technique  and  the  approximate  FFT  technique  go 
together  well.  The  octane  processing  and  variable  data  overlap  required 
b>  the  weighted-sum  approach  are  well  matched  to  the  approximate  FFT's 
characteristics.  The  benefits  of  the  combined  techniques  are 

1.  Computational  efficiency; 

2.  Variability  of  data  overlap;  anJ 

3.  Elimination  of  the  need  for  large  size  FFT's. 

There  are  some  difficulties  and  drawbacks  to  constant-Q  spectral 
analysis  by  means  of  the  combined  weighted-^um-approximate  FFT: 

1.  Different  data  rates  exist  within  and  between  octaves. 

2.  A picket-fence  effect  is  created  by  the  approximate  FFT  tech- 
nique. 
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3.  The  application  of  the  weighting  coefficients  is  complicated 
s lightly. 

The  existence  of  variable  data  rates  between  octaves  is  a funda- 
mental characteristic  of  constant-Q  processing.  For  example,  if  octave 
1 becomes  available  every  second,  octave  2 will  be  available  every  2 
seconds,  octave  3 every  4 seconds,  etc.  This  variable  data  rate  exists 
regardless  of  the  method  used  to  obtain  the  constant-Q  spectral  esti- 
mates. When  one  allows  for  variable  data  overlap  within  an  octave,  the 
data  rates  within  an  octave  also  become  variable.  There  is  no  easy  way 
around  this  complaint.  It  is  a problem  that  users  of  constant-Q  spect- 
tral  analysis  techniques  leam  to  accept.  An  example  of  the  variable 
data  rates  can  be  seen  in  the  next  section  of  this  report. 

The  picket-fence  effect  of  the  approximate  FFT  is  caused  by  the 
temporal  weighting  used  on  the  coarse  FFT.  This  effect  displays  itself 
when  one  looks  at  pure  tones  centered  on  vernier  FFT  outputs.  The  bin- 
centered  pure-tone  outputs  will  vary  in  amplitude  as  the  tone  varies 
from  bin-center  to  bin-center  within  the  vernier  FFT  of  a given  coarse 
spectral  bin.  This  variation  in  amplitude  is  known  and  can  be  corrected 
for  tones  centered  in  vernier  FFT  bins. 

Use  of  the  approximate  FFT  technique  as  an  octave  processor  with 
variable  overlap  of  the  temporal  FFT  outputs  complicates  the  application 
of  the  filter  weights.  Before  applying  a given  filter  weight,  one 
must  first  locate  the  uniform-resolution  bin  on  which  to  center  the  fil- 
ter weights.  In  the  case  of  the  approximate  FFT  technique  with  50  per- 
cent overlap  of  tne  temporal  FFT's,  one  must  first  determine  the  coarse 
spectral  bin  and  then  the  vernier  bin  on  which  to  center  the  filter 
weights . 

For  example,  suppose  we  had  to  center  a set  of  weights  about  f0  = 
256-1/16  Hz  and  had  performed  a length  8 vernier  FFT  on  coarse  spectral 
bins  255  and  256  Hz.  The  resulting  resolution  of  the  vernier  FFT  out- 
puts is  1/4  Hz  and  the  desired  fQ  lies  closest  to  the  zeroth  vernier  bin 
from  the  256  Hz  coarse  spectral  bin's  vernier  FFT  (see  figure  12).  Fig- 
ure 12  depicts  the  center  frequencies  of  vernier  FFT  outputs  rearranged 
by  considering  FFT  points  in  the  last  half  of  the  vernier  FFT  as  nega- 
tive frequencies.  There  are  redundant  spectral  estimates,  and  this  is 
the  reason  only  the  central  half  of  the  vernier  FFT  outputs  are  used. 

This  example  merely  illustrates  one  facet  of  the  implementation  of 
the  weighted-sum  constant-Q  spectral-analysis  technique.  Generally  one 
not  only  generates  sets  of  weights  representing  filters  within  an 
octave,  but  also  associates  with  each  set  a number  that  indicates  where 
the  filter  resides  within  an  octave. 
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THE  UNI VAC  1108  PROGRAM 


The  FORTRAN-based  Univac  1108  program  for  constant-Q  spectral 
analysis  by  means  of  the  approximate  FFT  was  originally  planned  to  be  a 
general-purpose  check-out  routine  for  the  technique.  Owing  to  inter- 
facing problems  between  existing  data  generation  and  display  routines, 
compromises  were  made  that  greatly  limit  the  routine's  flexibility. 

The  original  specifications  for  the  program  were 

1.  To  output  either  proportional  bandwidth  or  linear  spectral 
estimates; 

2.  To  perform  up  to  four  octaves  of  spectral  analysis;  and 

3.  To  divide  the  octaves  into  a maximum  of  four  sectors,  each  of 
which  could  have  a different  amount  of  overlap. 

What  currently  exists  is  a program  that 

1.  Outputs  proportional  bandwidth  spectral  estimates; 

2.  Performs  up  to  four  octaves  of  spectral  analysis;  and 

3.  Divides  each  octave  into  two  approximately  equal  halves,  where 
the  higher- frequency  half  has  75  percent  overlap  while  the  lower  has 

50  percent  overlap. 

Initial  outputs  of  the  program  are  very  encouraging  and  indicate 
that  the  merging  of  the  two  techniques  actually  works.  What  remains  to 
be  seen  is  whether  or  not  the  routine  is  as  efficient  computationally 
as  originally  envisioned. 

Figures  13  and  14  are  actual  normalized  outputs  from  the  program. 
Figure  14  is  a compressed  version  of  figure  13  and  displays  additional 
time  frames. 

In  the  figures,  the  horizontal  axis  represents  filter  numbers  (not 
center  frequency).  The  vertical  axis  represents  time  frames.  For  illu- 
stration purposes,  the  mapping  of  filter  number  to  frequency  and  time 
frame  to  seconds  is  irrelevant.  Table  1 gives  a breakdown  of  the  para- 
meters of  the  figures. 
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Figure  13. 
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Table  1. 


Octave  Setups 


Octave 
No . 

Sector 

No. 

Filters 

Percent 

Overlap 

Output  Period 
(Time  Frames/ 
Output ) 

1 

1 

580-750 

75 

1 

1 

2 

375-580 

50 

2 

7 

i. 

1 

215-375 

75 

2 

2 

2 

0-215 

50 

The 
fi  1 terin 

figures  illustrate 
g: 

the  followin 

g characteristics 

of  proportional 

1.  The  line  dynamics  are  compressed  because  of  the 
frequency  spacing  of  the  filter  center  frequencies;  and 

logarithmic 

2 # 

The  varying  overlap 

affects  the 

output  data  rates  of  the  sec- 

tors . 


The  effects  on  the  data  rate  can  be  seen  in  figure  13.  Note  that 
in  octave  1,  sector  1,  each  new  time  frame  has  a new  set  of  data.  In 
sector  2 of  octave  1,  the  new  data  set  occurs  every  other  time  frame.  By 
the  time  one  gets  to  octave  2,  sector  2,  the  new  data  set  occurs  every 
fourth  time  frame.  The  more  octaves  one  processes,  the  worse  this  data- 
output  condition  gets.  Another  interesting  effect  is  the  startup  tran- 
sient, which  also  gets  worse  as  one  processes  more  octaves.  This  startup 
transient  can  be  seen  in  the  number  of  time  frames  needed  to  obtain 
nonzero  data  in  a given  octave. 

The  compression  of  line  dynamics  can  be  seen  in  figure  14.  The 
lines  in  octave  1 are  harmonically  related  to  the  similarly  varying 
lines  in  octave  2 (e.g.,  the  line  around  filter  715  varies  twice  as 
much  in  frequency  as  the  corresponding  line  around  filter  340) . Note 
that  the  frequency  excursions,  in  terms  of  filter  numbers,  are  actually 
equal . 

In  order  to  generalize  the  program,  it  would  have  to  be  completely 
rewritten,  based  upon  the  basic  techniques  presented  in  earlier  sections. 
Things  such  as  data  structures  are  critical  to  the  efficiency  of  the 
technique,  and  perhaps  a more  efficient  method  of  implementing  the  over- 
lap (such  as  phase  shifting  the  filter  weights)  could  b^  devised. 
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SUMMARY  AND  RECOMMENDATIONS 


This  report  has  presented  a worthwhile  and  efficient  method  for 
obtaining  constant-Q  spectral  analysis,  along  with  some  of  the  pros  and 
cons  of  the  technique.  Also  represented  is  a powerful,  useful  applica- 
tion of  the  approximate  FFT  technique. 

What  lies  in  the  future  for  the  combination  of  these  techniques 
rests  in  the  hands  of  persons  interested  in  the  method.  The  author 
hopes  that  a valid  detailed  comparison  of  this  combination  of  techniques 
with  other  techniques,  with  respect  to  efficiency  and  accuracy,  will  he 
performed  eventually.  It  would  also  be  instructive  if  a comparison  of 
proportional-bandwidth  spectral  analysis  with  uniform-resolution  spec- 
tral analysis  were  performed,  to  see  if  the  gains  in  using  proportional 
bandwidth  are  worth  the  grief  of  struggling  with  variable  output  rates 
and  other  related  constant-Q  problems. 
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APPENDIXES 


SIDE  LOBE  CONTROL  FOR  FREQUENCY  DOMAIN 
SMOOTHING  IN  PROPORTIONAL  BANDWIDTH  PROCESSING 


Appendixes  A and  B originally  appeared  under  date  of  15  February 
1975  as  NUSC  Technical  Memorandum  No.  TC-4-75,  by  Albert  11.  Nuttall. 
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Appendix  A 

INTRODUCTION 

It  is  well  known  that  Hanning  temporal  weighting  on  a time  function  is 
equivalent  to  frequency  domain  smoothing  of  the  Discrete  Fourier  Transform 
(DFT)  coefficients  with  the  sequence  -1/4,  1/2,  -1/4.  However,  this  fixed 
smoothing  procedure  does  not  allow  any  variation  in  the  extent  (bandwidth) 
of  the  frequency  domain  smoothing  effected. 

If  the  time  duration  of  the  Hanning  temporal  weighting  were  contracted 
below  that  of  the  duration  of  time  over  which  the  DFT  were  conducted,  then  a 
broader  frequency  domain  smoothing  would  be  realized  However,  the  equivalent 
frequency  domain  procedure  would  be  smoothing  with  a long  and  more  complicated 
sequence  than  that  given  above.  It  is  the  purpose  of  this  memorandum  to 
investigate  the  possibility  of  using  frequency  domain  smoothing  via  short  sequences 
for  arbitrary  smoothing  bandwidths  and  to  point  out  any  limitations  of  the 
technique. 


TECHNICAL  CONSIDERATIONS 
Let  waveform  x(t)  have  voltage  density  spectrum* 

*ftr)  (i) 

Suppose  x(t)  is  sampled  at  times  where  a is  chosen  smaller  than  (2fH)  , 

and  is  the  highest  frequency  contained  in  X(f);  that  is,  no  aliasing  occurs. 

We  compute  the  DFT  of  the^available  samples: 

A*  = *(**)  0 - K * N-l.  ( z ) 

But  by  use  of  (1 ),  we  can  express 

A*  & WM Jf  «»j>(;2Tf»«)Xff) 
■JJFXW-f).  X(f)®-R(f)|  (3) 

*Time  functions  will  be  denoted  by  lower  case  symbols;  their  Fourier  transforms 
(frequency  functions)  will  be  denoted  by  capitals. 
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Equation  (3)  states  that  Fourier  coefficient  \ is  the  result  of  convolving  spectrum 
A(f)  with  window  R(f),  and  then  sampling  at  frequency  However,  R(f)  has 
bad  sidelobes  (v-13  dB  at  f=±^). 

In  an  attempt  to  reduce  sidelobes,  suppose  we  consider  the  sum  of  weighted 
Fourier  coefficients  (centered  about  some  frequency,  of  interest): 

tf-L 


C(Q * ir 

weights  fM^are  yet  to  be  specified.  This  is  the  frequency  domain  smoothing 
procedure  mentioned  in  the  Introduction.  Employing  (3),  (5)  becomes 

C(f.)  .£g^«>jdFX»)T?(fr-f) 


(s) 


ft 


where  window  function 


The  problem  now  is  to  choose  weightsfiU^uch  that  $ (f;  fo)  has  good  sidelobes 
versus  f,  in  the  neighborhood  of  fo.  If  we  accomplish  this,  (6)  indicates  that 
coefficient  C(fo)  is  the  result  of  looking  at  X(f)  through  a window  S(f;fo)  with 
good  sidelobe  behavior. 

To  get  at  the  choice  of  weights  consider  a real  and  even  weighting 
function  ^toof  duration  He.  Then  window  We  (f)  is  also  real  and  even;  see 
figure  1 . Define  a time-delayed  weighting 

wte)  * we(£--r),  (?) 
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I. 


Evan  WfrioVtlHtf 


9nJ  Window 


for  which  the  window 


WW)'  W,(f)  "fi-OrtT);  (?) 


see  figure  2.  The  magnitude  characteristics  and  sidelobes  of  W(f)  are  identical 

IvK)  «* — Na  — * 


to  those  of  We  (f). 


Figure  *•  Pel««jad  Wai^^Finj 


Now  consider  samples  of  a frequency-shifted  and  expanded  version  of  window 
W according  to 


v(^)  - 
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hi ow  if  u>*  )«c“t 

T * 

*here  I is  an  integer,  (11)  can  be  written  as 


Q*) 


-e^Ovf.A/ir).  (,3) 

That  is,  the  k -dependent  terms  are  real,  under  the  selection  (12)  for  delay  T 
The  choice  of  integer  I is  not  yet  determir^d. 


We  now  select  the  weights {B«ft)}to  be  real,  according  to  the  samples 


3,  If.)* 


, (») 


■tt 


others  se 


where 


Na 


<f  < 

< T®  hit 


by  this  choice  of  real  weights  for  (Bk(^)j. 

We  must  now  show  that  the  choice  (14)  results  in  a window  S(f;fo)  with  good 
sidelobe  properties.  Equation  (7)  yields  (using  realness  of  We) 

No*u  if  If,  K>  «"*  close*  SMtl  Uot 

Hwfo)l,  (n) 
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the  edge  terms  in  the  sum  in  (16)  are  very  small  compared  to  the  center  terms.  And 
if  window  We  (f)  has  a rapidly  decaying  sidelobe  behavior,  (16)  can  be 
approximated  by  the  infinite  sum 

*•(&-<) 

- i)Jf  V 

. -if  «,(»(-; x)[w(i^)L(f)} ® R*(- f),  (f *) 


where  impulse  train 


(It) 

The  Fourier  transform  of  S(f;fo)  is  then  available  from  (18)  as 

s (t)  f.)  2 § $ t)  S(f ) fo) 

- d ^-i  7rf0H4  *(«+) « r*t) 

s A i-rr^NA  x)  « *(* fe-kN*])  QXf>(i  M»&-  k (2o) 

where 

rlfrl-J*f  (i  ^ & (£-  »»*)  (2)) 

Notice  that  r(t)  is  non-zero  only  in  the  range  (Pj(A/-()a] 
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in 


Now  consider  the  plot  of  the  magnitude  of  the  sum  in  (20).  It  is  depicted 
figures  4A  and  4B  for  I even  and  I odd,  respectively;  we  have  employed 

F/jwr’c  4A.  X even  * * 

■- . A./|\  lAi.  A.,-.. 

f\  Ma  2Ma 


-til. 


Fiqure  4 9.  X <><}J  "l  ^ 

f\  A !A  A 

, / i v ! y / i v » 4 i a 

/N  NA  1.. 


-Na  6 Na  2AIA 

P»,guve  4.  of  5»w  in  (ao)  ( <Jn»w)n  "ftr  *<>  >). 

C-.miifr  2 «*»d  Gz>.  Th ft  Vu«oti*H  , _ 

^ • /or\\  . 1 1.1 «llnn  { ft  k)  In  fintirD  4 If  T It  PU«n  . t 


Tiaere  a «*»d  Gz).  The  function 

in  (20)  gates  out  the  portion  (0,A/a)  in  figure  4.  If  I is  even,  the 

resulting  function  sft)-f0)  has  large  discontinuities  at  its  edges  at  0 and  Ma, 
thereby  leading  to  a window  S( f ;£)  with  bad  sidelobes;  see  figure  4A.  But 

if  X is  odd  and  if  . » 

OC  > I,  (22^ 

figure  4B  shows  that  will  have  the  smooth  behavior  dictated  by  that 

shown  in  figure  2.  Hence  we  restrict  consideration  to  I odd  in  the  following. 

The  only  term  in  (20)  which  contributes  to  s(fc;£)  is  that  for  K*  -(r-iVi 
yields 


and  yi 


tt)  -Q  « a 

ex^iltr^^b-h  -TT-AJd) 

= (23) 

where  we  have  employed  (8)  and  (12).  Bui  since,  by  (22)  and  figure  1, 
*.(«&- 4?])  is  zero  outside  the  interval  (0j  AlA/  , we  con  rewrite 

(23),  using  (21)  and  (19),  as 


5 k.i)  * & S*i0#  H.Wt-50) 
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thereby  immediately  yielding 

S( (2r) 

The  origin  impulse  in  the  train  ^_l(0  yields  the  term 

We  VXy^i^f  M&),  (*0 

which  has  the  desired  magnitude  characteristics.  The  other  impulses  in 
lie  in  frequency  ranges  where  X»  is  zero  (by  proper  anti-aliasing  filtering 
prior  to  sampling  x(f)  every  4 seconds).  Hence 

S ff , •vxfl-i  W Na)  i^n) 

in  the  frequency  range  of  interest.  This  is  the  main  result  of  interest. 

An  alternative  expression  to  (6)  for  coefficient  C&)  is  available  from  (1) 

as 

C(0  = 

- y cH:  X (0  J <lf  S*(f ) i)  ft) 

But  from  (24)  and  figure  1, 

s*ll  jf.)  = A ^ It) < We(«fffc-  TT-])  «*/>(-!  *"6 [t- ^]) . 

The  complex  exponential  in  (29)  " beats  down"  the  spectrum  of  xfc)  by  4; 
the  weighting  function  w*  gates  the  portion  of  x|fc)  indicated  in  figure  5;  and 
the  impulse  train  samples  this  gated  function.  The  integral  of  the  product 
of  X OwA  S is  then  taken  according  to  (28).  Figure  5 and  (28)  illustrate 
that  the  edges  of  waveform  x/t)  (near  OaW*>)  will  receive  little  or  no  weight, 
especially  if  V > I ; this  feature  may  require  overlapped  processing  of  *(£•). 


(7%) 

fa) 
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Notice  that  window  SIM)  in  (27)  i$  independent  of  X (so  long  os 
X is  odd).  Hence  we  can  choose  X * I in  (14)  with  no  loss  in  generality. 

Notice  also  that  the  sidelobe  characteristics  of  are  exactly  those  of 

V«.  irrespective  of  shift  -fo  and  expansion  •<  , provided  that  and  that 

(17)  is  satisfied.  Thus  we  have  merely  to  ensure  that  we  select  a good  window 
in  the  first  place  when  we  start  computing  weights  via  (14),  in  order  to 
guarantee  a good  window  of  course  "Til  is  also  necessary. 

There  is  an  upper  limit  on  e< ; that  ^arbitrarily  large  expansion  factors  in 
frequency  are  not  allowable.  Since  K,  and  |fa  are  bounded  as  in  (14) 
and  (15),  there  will  come  a value  of  < beyond  which  the  skirts  of  ^ are  not 
being  sampled  adequately.  The  approximation  in  (18)  is  then  inadequate  and 
the  succeeding  results  inapplicable.  For  fixed  If,  and  K2,  very  stringent  upper 
limits  on  «(  apply.  For  example,  it  is  found  that  for  Ka~)C,  * 6 (7  terms  in 

(14)),  values  of  •(  in  the  range  , \ 

1 < *C  S 2 (3° ) 

are  allowed  for  sidelobes  in  the  -30  dB  to  -40dB  range.  To  go  beyond  «f  = 2- 
will  require  more  weights  than  7.  Also  additional  overlapped  processing  will 
be  required;  recall  figure  5. 

The  restriction  to  finite  K,  and  K,  in  (1 4)— (1 6)  (rather  than  the  infinite 
assumption  of  (18))  makes  the  sidelobes  of  3fF)£)  depend  on  ^ , rather  than 
the  simple  dependence  in  (27).  However,  these  sidelobes  can  be  controlled  by 
a sufficiently  large  value  of  ICj-ki . 

Fih»ll>,  (2t)  , (c)  cam  !>«  e*prtjse^ 
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C(£)  * « A,  3,(1) 

^ X(f)  W,  (ti)  (a/) 

The  phase  factors  ir  (4)  and  (27)  are  of  no  concern  as  regards  their  effects  in 
(3)  and  (6)  respectively,  >ince  they  can  be  absorbed  by  a redefinition  of  the 
time  origin  of  xft)  in  both  cases. 


DISCUSSION 

The  desired  sideloh**  behavior  and  variable  bandwidth  expansion  factor 
can  be  realized  through  frequency  domain  smoothing,  provided  a sufficient 
number  of  coefficients  are  utilized  in  the  frequency  domain.  However,  the 
larger  bandwidth  expansion  factors  will  require  the  use  of  more  coefficients. 
Additionally,  the  equivalent  time  domain  operation  of  multiplication  will 
squelch  the  edges  of  the  data  and  may  thereby  require  employment  of  overlapped 
processing.  These  limitations  must  be  kept  in  mind  when  this  technique  is 
utilized. 
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Appendix  B 


PROGRAM 


Letf  «*;  + > 05  rs  10. 

(Let  r*<0  in  order  to  look  to  left  of  ).  (l)t 

S9iQ' 


Let 

and 

Then 


k,*0,  - 2M. 

N-l 

2 

k=o 


r~t  M-K-  m,\ 
NA  ) 


C*l) 

(*$ 

(A3) 

(Av) 


But  from  (4), 


-W/JC 


w * W ' "rv " -w-‘ y N 9i„  (^) 

Therefore 

5(hi)  = * £(«■)  (t<) 

|C  20 


Also 


and 


Therefore 


3„(f.)  *(-'?  H(^-)  «t)  * Q», , °iKs  ^ 

(at) 

c,tMi  =(-') (r?1- *),  - M - K 1 M ■ (p^r  (a*) 

Grid  = 
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For  the  Kaiser-Bessel  window, 


m')  - 

tjff-  tlrbl ) 


C = (- i)  Fws(>) 


- h <k<  M 


where 


5'ihV»  (\/ £x~x'>) 


5 • 


For  the  Hanning  window, 


FNS(*)  = 


5 »n  x 

*0-(SO 


A basic  program  for  the  HP9830A  follows. 


(<  lo) 

(am) 
(a  ix) 

M 
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420  IiEF  FNR'Rl  ' 

430  IF  R1--0  THEN  480 
440  R2=P I *R 1 
450  R3=R2/N 

460  R=C08  ( ( N- 1 ) *R3  > * S I N < R2  1 *.  N+s  1 N F'  •.  •. 

470  RETURN  R 
430  R- 1 
490  RETURN  R 
500  DEF  F N I < R 1 > 

510  IF  R 1 — 0 THEN  560 
520  R2=PI *R1 
530  R3=R2  H 

54U  I = 3 1 N < N 1 ) * R 3 > * o 1 H ( F'  2 1 N *■  . I (1  '•  R 1 1 

550  RETURN  I 
560  1=0 
570  RETURN  1 
530  DEF  FHSOO 

590  REM  K RISER  BE 3 3 EL s G I HH 0 1 11  NHFRE 

600  REM  HRHN I NG : S I N X > ■■  v*  < i-( : : p|  .t  -■  . 

610  S 1 =B 1 2-X  t'2 

620  IF  SI  •■=  O THEN  660 

630  S1=SQR'-S1) 

640  S=S I N 1 SI  ' SI 

650  RETURN  S 

660  IF  SI  0 THEN  690 

670  3=1 

630  P E T U R N 

690  3 1 =SQR < S 1 • 

700  S2=EXP *;  SI  1 
7 1 0 S = (.  S 2 ~ 1 3 2 1 • 2 * S 1 1 
720  RETURN  S 
730  END 
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